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I.   INTRODUCTION 

Twenty  years  ago  discussions  of  the  relative  merits  of  multistep 
versus  RK  methods  made  frequent  reference  to  the  "difficulty"  of  starting 
a  multistep  method;  that  is,  the  cost  of  obtaining  the  additional  k  -  1 
values  needed  to  start  a  k-step  method  when  just  the  initial  value  of  the 
problem  is  available.   Some  early  programs  for  constant  step,  constant 
order,  multistep  methods  used  several  RK  steps  to  start.   For  example, 
a  four-step  Adams  formula  can  be  started  with  three  fourth  order  RK 
steps  to  obtain  y!,  i  =  0,  1,  2,  3  and  y  from  y  .   No  error  control  was 
provided  in  such  a  program,  but  at  that  time  the  error  control  in  most 
programs  consisted  of  an  external  adjustment  of  the  step  by  the  user. 
Nordsieck  [12]  designed  one  of  the  first  multistep  methods  with  automatic 
starting,  but  the  introduction  of  variable-order  methods  finally  "solved" 
the  problem  because  these  methods  can  start  at  whichever  order  corresponds 
to  a  one-step  method — usually  first  or  second  order.   See,  for  example, 
Gear  [6],  Hindmarsh  [9],  Krogh  [11],  Shampine  and  Gordon  [13].   (It  is 
interesting  to  observe  that  Nordsieck' s  scheme,  which  consisted  of 
integrating  forwards  and  backwards  over  several  steps  several  times, 
amounts  to  a  variable-order  scheme:   because  Adam's  methods  are  used  and 
because  of  the  representation  used,  it  can  be  shown  that  the  order  is 
increased  at  each  step  in  the  Nordsieck  starter.) 

Starting  at  first  or  second  order  and  relying  on  the  order 
control  to  increase  the  order  to  an  appropriate  value  is  programmatically 
simple  because  the  mechanism  is  already  present  in  the  program.   It  also 
has  the  advantage  of  providing  local  error  control,  which  was  not  really 
present  in  the  Nordsieck  scheme.   However,  it  is  not  very  efficient  for 
any  but  low  accuracies  as  very  small  stepsizes  must  be  taken  at  low  orders 


to  maintain  accuracy.   This  inefficiency  is  often  offset  by  the  improved 
efficiency  of  high  order  multistep  methods  for  many  classes  of  problems 
(see  Hull  et  al.  [10]).   However,  there  are  large  classes  of  problems, 
arising,  for  example,  in  simulation,  in  which  there  are  frequent 
discontinuities  in  first  or  higher  derivatives.   Since  multistep  methods 
are  based  on  polynomial  interpolations  over  an  interval  of  several  steps, 
they  break  down  when  the  interval  includes  a  discontinuity  in  a  derivative 
of  degree  lower  than  the  polynomial.   Hence,  a  restart  is  necessary. 
One-step  methods  such  as  RK  do  not  have  this  difficulty  as  long  as  a  mesh 
point  falls  exactly  on  each  discontinuity,  and,  therefore,  RK  methods  are 
usually  superior  for  these  "non-smooth"  problems. 

This  paper  presents  a  RK-like  technique  for  starting  (or 
restarting)  any  of  the  currently  popular  automatic  multistep  ODE  integrators 
at  a  high  order.   (An  order  four  starter  will  be  recommended  for  general 
use.)   The  technique  uses  fewer  function  evaluations  than  are  used  if  a 
conventional  RK  method  is  applied  repeatedly,  and  provides  some  degree  of 
error  control,  as  well  as  data  for  estimating  the  step  size  to  be  used  in 
the  first  step  of  the  multistep  method. 

Section  2  presents  the  idea  in  a  non-RK  framework  to  establish 
the  existence  of  methods  of  the  desired  type  and  to  derive  an  upper  bound 
on  the  number  of  function  evaluations  needed  in  an  explicit  RK  starter 
(1  +  p(p  -  l)/2  for  a  p-th  order  method).   Section  3  examines  the  lower 
bound  on  the  number  of  function  evaluations  in  explicit  RK  starters  for 
p  <^  4,  and  suggests  a  particular  fourth  order  method  that  uses  six 
function  evaluations  (the  lower  bound).   Section  4  examines  implicit  RK 
starters  which  might  be  useful  for  some  stiff  problems  and  proves  that 
the  minimum  number  of  stages  for  a  p-th  order  starter  is  p,  and  that  p-stage 
methods  with  desirable  stability  properties  exist.   Section  5  considers 


programming  considerations  and  error  control,  while  section  6  discusses 
the  detection  of  discontinuities. 


2.   STARTING  BY  EXTRAPOLATION 

If  a  k-step  multistep  method  is  to  be  started  at  p-th  order,  it 
is  desirable  to  generate  k  -  1  additional  starting  values  which  are  p-th 
order  accurate.*   Since  k  is  p  -  1  or  p  in  most  methods  in  use  (Adams 
or  backward  differentiation  formulas),  we  will  consider  the  case  k  =  p. 
If  we  can  generate  hqy    ,~, ,  q  =  0,  1,...,  p,  with  error  no  greater  than 
0(h   ) ,  we  can  easily  compute  the  starting  values  needed  for  a  code  which 
stores  either  y   (approximations  to  y(t.)),  y! ,  or  h^y^q  /q!  (the 
Nordsieck  vector).   One  way  of  doing  this  is  by  extrapolation  as  shown 
below.   (It  should  be  emphasized  that  the  following  is  an  existence 

argument.   It  is  not  the  way  to  compute  the  required  values!) 

i 
Let  y.,  i  =  0,...,  £  be  the  results  obtained  from  the  application 

of  the  forward  Euler  method  for  £  steps  to  the  differential  equation 

y'  =  f(y,  t) ,  y(0)  =  y_,  using  step  size  h  =  H/&.   We  will  assume  that 

f(y,  t)  has  as  many  derivatives  as  needed  (for  example,  bounded  fourth 

partial  derivatives  will  be  needed  for  a  fourth  order  method).   Then  it 

is  known  that 

(2.1)  yf=y(ih)+  I     hqe  (ih)  +  0(hP+1) 

q=l    q 

(see  Stetter  [14].   Also,  there  exist  sets  of  numerical  differentiation 
formulas  for  k  +  1  equally  spaced  points  of  the  form 

(2.2)  Z  d.,z(ih)  =  hkz(k)(kh/2)  +   Z   c  hSz(s)(kh/2)  +  0(hP+1) 
i=0  ±tC  s=k+l  SlC 

(In  fact,  c  ,  =  0  when  s  -  k  is  even.) 


*  If  restarting  occurs  a  fixed  number  of  times  in  a  fixed  interval, 
errors  of  order  hP  in  the  starter  do  not  cause  a  total  error  of  more  than 
hP  in  the  global  solution,  but  a  practical  algorithm  may  restart  more  and 
more  frequently  as  the  error  tolerance  is  decreased,  so  we  will  try  and 
keep  the  restarting  error  to  0(hP+l). 


Applying  eq.  (2.2)  to  eq.  (2.1)  with  k  =  p,  h  =  h  =  H/p  we  get 


(2.3) 


S  d.  y?  =  hPy(p)(H/2)  +0(hP+1) 
1-0  ip  l         p  p 


Since  the  left-hand  side  of  eq.  (2.3)  can  be  evaluated,  we  can  compute  an 
order  p  accurate  approximation  to  H  y    .   From  now  on,  we  will  omit  the 
term  0(h   )  to  ease  the  eye  strain. 

Now  apply  (2.2)  to  (2.1)  with  k  =  p  -  1  and  h  =  h  1  and  h  . 


We  get 
(2.4a) 


PZ  d.     y?"1  =  hP"^  y(P_i;(H/2)  +  hP  ,  ejp-i;(H/2)  +  c„  1 
1=0 


p-1  _  p-1  „(p-l) 
'i  p-1  yi      p-1 


(P-D 
p-1  "I 


(P) 


PP  P+l 


(2.4b) 


p-1 


Ed.    1  yP  =  h15"1  y(p_1)(H/2  -  h  /2)  +  hP  e1(p"1)  (H/2  -  h  /2) 

i=Q   X  P-1   !      P  P  P    1  P 


+  c   hP  y(p)(H/2  -  h  II) 
PP  P  P 

Note  that  h?-1   y(p-1)(H/2  -  h  ID   -  h^1  y(p"1)(H/2)  -  \   hP  y(p)(H/2) 
P  P        P  2   p 

and  hP  e1(p_1)(H/2  -  h  /2)  =  hP  en(p_1)  (H/2)  (neglecting  terms  of  order  hP+1). 
pi  p       pi 

Furthermore,  we  know  h  y    (H/2)  from  eq.  (2.3)  to  order  p,  hence  eqs.  (2.4) 
can  be  rewritten  as 

hP-l  y(P-1)(H/2)  +  hP  ,  en(p_1)(H/2)  =  known  values 
p-1  p-1   1 


and 

hP-l  y(P-1)(H/2)  +  hP  e1(p_1)(H/2)  =  known  values. 
P  P   1 

These  can  be  solved  for  HP_1  y(p_1)(H/2)  and  HP  e^P_1) (H/2)  to  order  hP 

since  h  ^  h 

P    P-1 

The  process  continues  in  a  similar  fashion.   Suppose  we  have 
already  computed  HS  y(s)  (H/2)  for  s  =  k  +  1, .  .  .  ,  p,  and  H^"*"3  efs)(H/2)  for 
j=l,...,p-k-l,  s=k+l,...,p-j.   We  now  apply  (2.2)  to  (2.1) 


with  h  =  h.  ,  h.  ,,,...  ,  and  h  .   We  get 

K.     K.+1  p 

(2.5a)     I  d.k  yk  =  hk  y(k)(H/2)  +  h£+1  e<k)(H/2)  +.  .  .+  h£  e(^(H/2) 
i=0  ^ 

(2-5b)    ,*n  dik  yk+1  ■  hu+i  *(k)(H/2  -  hk+i/2)  +  C  eik>(H/2  -  \+i/2)  +- 

1=0 


k 

E 
i=0 


(2.59)     E  d.,  y?  =  hk  y(k)(H/2  -  (p-k)h  /2)  +  hk+1  en(k)(H/2  -  (p-k)h  /2)  +. 
,  rt  ik  i    p  p      pi  p 


(k)       (k) 
As  before,  the  values  of  y    and  e    at  H/2  -  mh  can  be  shifted  to  values 

at  H/2  by  means  of  Taylors  series  using  the  already  computed  higher 

derivatives  of  y  and  e  .   After  this  transformation,  eqs.  (2.5)  are  a  set 

q 

(k) 
of  p  +  1  -  k  linear  equations  in  the  P  +  1  -  k  unknowns  y    (H/2), 

(k)  (k) 

e1   (H/2),...,  e   ,  (H/2)  whose  coefficient  matrix  is  a  modification  of  the 
1  p-k 

Vandermonde  matrix,  and  is  nonsingular  since  the  h  's  are  all  different. 
Thus,  this  process  yields  order  p  correct  values  to  H4  y    (H/2)  for 
q  =  1,  2,...,  p.   The  "cost"  of  this  in  function  evaluations  can  be  seen  to 
be  1  +  p(p-l)/2  because  f(yn»  0)  is  evaluated  once  and  used  in  all  p 
integrations.   An  example  is  presented  at  the  start  of  the  next  section. 


3.   EXPLICIT  RUNGE-KUTTA  METHODS 

The  extrapolation  technique  of  the  previous  section  can  be  viewed 
as  a  Runge-Kutta  method.  We  will  first  illustrate  this  for  the  case  p  =  3. 
Let  h  =  H/3  and  consider  an  autonomous  system.   We  have 

r  3 


(3.1)       ( 


yx   =y0  +  hf(y0)  =y0  +  kQ 


yJ-yJ  +  hfCy?)  -  y0  +  %  +  \ 


jl-732  +  mj32)  =y0  +  ko  +  ki  +  k2 


yl  =  yo+Tf(V  =  yo  +  2  ko 


where  kQ  =  hf (yQ) 


3 
where  k  =  hf (y  ) 


where  k  =  hf(y_) 


2    23h£/2v       ■  3  .   ,  3  , 

y2  -  yj+T^)  =  y0  +  2  ko+  I  k3 


k. 


yj  -  y0  +  3hf (yQ)  =  yQ  + 


3k. 


where  k  =  hf (y  ) 


We  use  the  difference  formula  given  in  (3.2)  below  where  all  derivatives 

4 
are  evaluated  at  t  =  3h/2.   0(h)  terms  are  neglected. 

"3    3     3   _  3        3   (3) 

D3  -  y3  -  3y2  +  3Yl  -  yQ  =  h  y 

i  n3    3    3    3        ,,2   (2)     3   (2) 

2  =  y3  "  y2  "  yl   y0       y         el 


(3.2) 


:  n2  2     2        3h.2   (2)    ,3hv3   (2) 

J  D2  =  y2  "  2yl  +  y0  =  (T}   y    +  (T)   el 

t 

n3  3    3   .  (1)   .2   (1)  ,  .3   (1)   h3   (3) 

1  y2  "  yl  =  y         el        e2     24  Y 


(1)  ,  9^2  _(1)  ^27^3  _(1)  ^9^3  „(3) 

1 


Dl  =  y2  "  y0  =  3hy    +  2h  el   +  Th  e2   +  «   Y 


nl    1        iT,  (1)    Q,2   (1)    ._,  3   (1)  ,  27  ,3   (3) 

^-  1  =  Yl  "  y0  =   y  el  e2     24    y 


(Symmetric  differences  have  been  used  above  to  reduce  arithmetic.)   From 
these  we  find  that 


namely 


(3.4) 


'V  y(3>   -  D]  +  0(hS 


(3.3)        I   h2  y<«   -  |  D3  -  |  »\  +  o(h4) 

9n3       4  ^2   J    1  _1        9  „3 


2DriDl  +  6Dl+lV°(h) 


Fairly  clearly,  we  can  view  this  as  a  Runge-Kutta  process, 


JL 

f  k.  =  hf(yn  +  Z   3..  k.   ), 
i   i       0    ml     ij   j-1 


i  =  0,.. . ,  q-1 


LhS  y(s)  =   Z  k-_!T.s+  0(hP+1),     s  =  1,...,  p 
where,  in  the  example  above,  q  is  four,  the  matrix  of  3  coefficients  is 


B  = 


0 

0 

0 

0 

1 

0 

0 

0 

1 

1 

0 

0 

3 

0 

0 

0 

and  the  matrix  of  y  coefficients  is 


3/2 


"-3/8   -1/6    1 
9/4     0     -2 
9/8    3/2    1 
_  -2    -4/3    0_ 

This  gives  estimates  for  the  derivatives  at  H/2  =  3h/2.   Estimates  of  the 
derivatives  and/or  function  values  to  the  same  order  of  accuracy  can  be 
obtained  at  any  point  within  a  constant  multiple  of  the  interval  H.   We 
used  the  midpoint  above  to  take  advantage  of  symmetry,  but  for  convenience 
in  the  discussion  below  we  will  consider  the  problem  of  computing  the 


derivatives  at  the  origin.   In  that  case  the  matrix  above  becomes 


1   -5/3 

0     3 
r°=   0    0    1 

0   -1/3    0 
Obviously,  the  estimate  for  hy  (0)  is  simply  k  . 

In  the  example  above  it  took  A  function  evaluations  for 
third  order.   Using  the  technique  of  section  2,  it  would  take  7  function 
evaluations  for  fourth  order  and  11  for  fifth.   It  is  well  known  that 
Runge-Kutta  methods  of  third,  fourth,  and  fifth  orders  are  possible  with 
3,  A,  and  6  function  evaluations,  respectively.   This  naturally  raises  the 
question:   Can  the  first  p  derivatives  be  found  to  accuracy  0(h    )  with 
fewer  than  1  +  p(p-l)/2  function  evaluations?  Note  that  this  is  not  the 
same  problem  as  finding  embedded  RK  methods  of  several  orders  such  as 
the  RK  Fehlberg  method.   See  Bettis  [1].   For  p  <^  3  the  answer  is  no, 
but  for  p  =  A  it  is  possible  in  six  but  no  less.   The  cases  p  =  1  and 
p  =  2  are  trivial.   The  problem  is  examined  in  detail  for  p  =  3  below 
and  for  p  =  A  in  the  appendix. 


3.1.   Non-Existence  of  3-Stage  Method  of  Third  Order 

In  the  discussion  below,  all  variables  are  evaluated  at 
t  =  0,  y  =  y  ,  unless  specified  otherwise.   Thus  we  can  write 
hy'  =  hf  =  k  .   For  a  system  of  m  equations  in  y  ,  y  ,...,  y 


(3.5) 


h2y"  =  h2   I^f1 


i=l  dyX 


10 


We  will  write  the  right-hand  side  as  h  f-,f.   Similarly 


(3.6) 


3  m  =  h3    ™    (jh_r   l^^-^Kfh 


h  y 


i,j=l   8yX8yJ 


3y  9yJ 


=  h3(f"2f2  +  f2f)     (definition  of  f£) 


Writing  ot.  =  £   3..,  the  k.  in  eqs.  (3.4)  can  be  expanded  to 

k0  =  hf 


ot 
k  =  hf  +  a  h2f1f  +  -y  h3f2f2  +  0(h4) 


ct 
k2  =  hf  +  a2  g2f1f  +  -y  h3f2f2  +  ax  322  h3f2f  +  0(h4) 

Therefore,  the  second  of  eqs.  (3.4)  can  be  written  as 


(3.7) 


1 

1 

1 

ni 

0 

d 

0 

al 

a2 

r  = 

0 

i 

0 

0 

a2/2 

a2/ 2 

0 

0 

l 

0 

0 

aie22_ 

0 

0 

l 

where  T   is  the  matrix  [y.  ].   The  rows  of  this  equation  correspond  to 

is 

2       2 
the  terms  f,  f,f,  f9f  ,  and  f,f>  respectively.   Clearly,  the  first  column 

T 
of  T  is  [1,  0,  0]   and  the  first  row  is  such  that  the  sum  of  all  rows  is 

[1,  0,  0],   Thus,  the  first  row  and  column  can  be  calculated  if  values 

can  be  found  to  satisfy  the  remaining  six  equations  in  eq.  (3.7). 

Therefore,  we  will  drop  the  first  row  and  column  in  (3.7)  to  get 


11 


(3.8) 


a 


a. 


a2 

2 
a2 

r  = 

ai322_ 

where  T  is  a  2  x  2  matrix.   If  this  is  to  have  a  solution,  a  S   ^  0. 
The  equation  in  the  (3,  1)  position  of  eq.  (3.8)  implies  Y91  =  0. 
Then  the  equation  in  the  (2,  1)  position  implies  Y, ,  =  0  so  that  is 
impossible  to  satisfy  the  equation  in  the  (1,  1)  position.   Hence,  a 
3-stage  third  order  method  does  not  exist.   However,  we  have  demonstrated 
the  existence  of  a  4-stage  method. 


3.2.   Fourth  Order  Methods 

Equations  (3.4)  must  be  satisfied  as  identities  in  each  of  the 
elementary  differentials  of  f  of  orders  up  to  p.   This  leads  to  a  large 
system  of  nonlinear  equations  for  large  p.   The  number  of  elementary 
differentials  for  p  <  5  are  given  in  Table  3.1. 


TABLE  3.1.   Elementary  Differentials  of  f 


Order    Number 


Representation 


fif 

2  2 

V  •  £if 

3  2       2   3 
f3f  ,  f2fxf  ,  flf2f  ,  fxf 

f4f4'  f3flf3'  flf3f3'  f2f3'  f2flf2> 
2         2   2   2   4  4 

f2(f1f)z,  fjfjfjf  ,  qf2f  ,  fjf 


The  elementary  differentials  are  represented  in  a  prefix  operator  notation: 
for  example,  f.  is  a  ternary  operator  given  by 


12 


m   m   m 
f~  abc  = 


3  f     ij   k 
-  a  be 


•  i  •  i  i  t  3y.3y.3y, 

i=l  j=l  k=l  'i  J2      k 
where  m  is  the  dimension  of  the  system  and  a,  b,  and  c  are  the  three 
vector  operands  of  f~.   For  an  alternate  notation  due  to  Butcher,  see 
Stetter  [15,  p.  Ill  ff.].   There  are  8  elementary  differentials  of  order 
up  to  four.   Thus,  the  second  of  eqs.  (3.4)  can  be  written  as 


a  r  = 


1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

3 

0 

0 

0 

1 

_o   o   o   : 

2   2      3       2       2 
where  the  rows  correspond  to  f,  f-,f>  f^f  >  f-,f>  f-jf  >  fo^i^  >  f-i^of  »  an(^ 


■2~r 


12 


f,f,  respectively.   For  a  q-stage  method  A  is  an  8  by  q  matrix  whose 
entries  are  completely  determined  by  the  $..,  and  T  is  a  q  by  4  matrix. 
As  before,  the  first  row  and  column  of  all  matrices  can  be  discarded 

because  the  first  row  of  A  is  [1,  1,...,  1]  and  the  first  column  is 

T 
[1,  0,...,  0]  .   This  leaves  a  system  of  21  nonlinear  equations  in  3(q-l) 

+  q(q-l)/2  unknowns.   Counting  unknowns  is  of  no  direct  value  in 

determining  the  existence  of  solutions,  although  it  can  be  a  guide  to  the 

prospects.   Although  for  q  =  5  there  are  21  equations  in  22  unknowns,  it 

is  shown  in  the  appendix  to  the  report  [8]  that  no  solution  exists. 

However,  for  q  =  6,  when  there  are  21  equations  in  30  unknowns,  there 

exists  a  9  parameter  family  of  solutions.   This  is  given  in  the  appendix. 

A  particular  case  of  this  is 
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B  = 


0 

0 

0 

1 

0 

0 

0 

2 

0 

0 

o- 

0 

0 

0 

0 

0 

0 

2 

0 

/3 

2_ 

,    r  = 


1 

-5/6 

4/9 

-1/9 

0 

0 

0 

0 

0 

1/2 

-4/9 

1/9 

0 

7/3 

-19/9 

7/9 

0 

-3 

10/3 

-4/3 

0 

1 

-11/9 

5/9 

3/4   0   9/4 

1/2   1   1/2 
1/12   2   1/4   2/3 

Naturally,  one  wonders  about  higher-order  techniques.   For 
example,  since  there  are  17  elementary  differentials  of  orders  up  to 
5,  we  are  led  to  a  system  of  16  by  4  equations  (after  discarding  the 
first  row  and  column)  in  (q+8)(q-l)/2  unknowns.   Nine  is  the  first  value 
of  q  for  which  the  number  of  unknowns  exceeds  the  number  of  equations, 
but  it  is  not  known  if  a  solution  exists  for  q  =  9.   Neither  does  it 
seem  practically  important  since  the  increase  in  the  amount  of  work 
does  not  seem  to  justify  the  small  increase  in  order. 
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4.   IMPLICIT  RUNGE-KUTTA  METHODS  AND  STIFF  EQUATIONS 

It  is  well  known  that  explicit  methods  do  not  function  well  for 
stiff  equations.   However,  normally  immediately  after  a  discontinuity 
a  differential  equation  is  not  stiff  because  the  solution  is  in  a  transient 
region  where  the  step  must  be  small  to  control  stability.   This  can  be 
seen  by  considering  the  equation 

(4.1)  y1  =  A(y  -  F(t))  +  g(t)  where  g  =  F'  almost  everywhere  and  X  <<  0. 
Its  solution  is 

(4.2)  y  =  ceAt  +  F(t) 

After  a  small  time,  the  terms  ce  '  is  small  compared  to  F,  so  the  behavior 
of  F  dominates  the  solution.   However,  if  there  is  a  discontinuity  in  the 
right-hand  side  of  (4.1)  at  some  £,  and  this  discontinuity  is  not  consistent 
with  g  =  F' ,  then  the  solution  (4.2)  will  have  a  discontinuity  in  c  and  F, 
probably  making  the  first  term  significant  for  another  short  time.   If  this 
is  the  case,  the  explicit  methods  discussed  earlier  can  be  used.   However, 
in  some  cases  the  discontinuity  may  not  stimulate  the  rapidly  decaying 
components  sufficiently  so  that  in  principle  a  step  size  large  compared 
to  1/X  can  be  used.   If  the  explicit  methods  Of  the  previous  sections  are 
applied  to  (4.1),  we  will  obtain  approximations  to  the  scaled  derivatives 
of  F  plus  polynomials  in  hX.   For  example,  if  the  6  stage  method  of  order 
4  is  used,  the  approximation  to  the  p-th  derivative  of  y  will  be 

hP  y<P>  =  hP  F<P>  +  c(hA)P  +  0(h5)  +  cefhX)5  +  ce  JhX)6 

lb  pb 

where  the  0(h)  term  depends  on  derivatives  of  F  only.   If  c  is  small, 

we  would  like  the  derivatives  of  F  to  dominate  this  expression.   However, 

if  hX  is  large,  the  (hX)   and  (hX)   terms  may  well  dominate.   For  these 

problems  it  is  worth  considering  implicit  RK  methods.   It  is  known 
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(Butcher  [2])  that  q-stage  implicit  RK  methods  can  achieve  order  2q. 
It  is  also  known  (Ehle  [4])  that  such  methods  are  A-stable.   In  fact, 
we  are  not  concerned  with  A-stability  or  variants,  such  as  stiff- 
stability,  here  because  only  a  single  RK  step  is  to  be  used.   What  we 
would  like  is  approximations  to  the  first  p  derivatives  of  y  to  accuracy 
h  with  the  property  that  they  are  not  large  for  hA  in  the  negative 
half  plane.   However,  note  that  we  would  like  these  approximations  to  be 
nonzero  as  hA  ■*■  °°  so  that  if  there  is,  in  fact,  a  significant  ce 
component  present  in  the  solution,  it  can  be  detected  and  the  stepsize 
used  in  the  starter  reduced. 

It  is  known  (Gear  [7])  that  a  q-stage  implicit  RK  scheme 
applied  to  y'  =  Ay  leads  to  the  calculation  of  terms  of  the  form 

(4.3)  R(hA)y=^Mly 

where  the  degree  of  the  polynomials  P  and  Q  cannot  exceed  q.   Since  we 

wish  to  get  an  order  p  approximation  to  h  y    ,  we  must  be  able  to  compute 

a  term  of  the  form  (hA)P  +  0(hA)P   .   This  implies  that  q  _>  p.   In  fact, 

we  will  demonstrate  that  it  is  possible  to  get  all  the  desired  terms  with 

q  =  p,  in  which  case  we  will  be  calculating  the  rationals 

2  p 

y  +  a  y  +.. .+  a  .  yF 

Rx(y)  = l ^— 

1  +  a-  y  +. .  .+  a  y 
1         P 

2      3  p 

y  +  ay  +...+  a    yK 

R2(y)  =  1 ^— 

(4.4)  1  +  an  y  +.  .  .+  a  yP 

1  P 


R  (y)  = 


1  +  a,  y  +.. .+  a  yP 
I         P 
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Note  that  for  these  approximations  to  be  bounded  when  U  is  in  the  negative 

half  plane,  and  to  be  nonzero  when  y  ■+  °°,  all  of  the  a  must  be  nonzero 

P 

and  all  roots  of  Q(u)  =  1  +  a  u  +. . .+  a  u  must  be  in  the  positive  half 

plane.   (It  is  probably  sufficient  that  a  ^  0  plus  the  root  condition  so 

that  at  least  the  approximation  to  h  y    is  nonzero  when  the  stiff 

components  are  significant.    In  fact,  a  should  not  be  too  large,  so  that 

a  large  value  of  h  y    can  be  used  as  a  warning  to  reduce  the  step  and 

try  again.) 

The  existence  of  p-stage  order  p  methods  for  starting  can  be 

demonstrated  without  reference  to  elementary  differentials — fortunately! 

Let  us  assume  that  £.  are  the  exact  values  of  hf(y(t.  +  ot.h)).   Consider 

i  0    l 

P 
k.  =  hf(y(t  )  +  Z  $   k.)    i  =  1,...,  p. 
1        u    j=i  ^   J 

If  the  a.  are  all  different,  the  interpolation   coefficients  3..  can  be 
chosen  so  that 


y(t0)  +  I  S   k 

-,=1     J    -J 

is  an  approximation  to  y(t^  +  a.h)  with  an  error  of  0(h).   Then,  k.  is 

Ox  i 

^  p+1 

an  approximation  to  k.  with  error  of  0(h   ).   Therefore,  solution  of 

P 
(4.5)        k.  =  hf(y(t  )  +  Z  3..  k.)    1=1,...,  p. 

will  lead  to  approximations  to  k.  with  errors  of  0(h   ).   With  these 

approximations  to  hy'  at  p  different  points  we  can  obtain  approximations 

s   (s)  P+1 

toh  y,  s=l,...,p  with  errors  of  0(h    ).   (In  fact,  determination 

of  the  interpolation  coefficients  3..  leaves  one  degree  of  freedom  in  each 
set,  which,  together  with  the  free  choice  of  the  a.'s  leaves  2q  degrees 
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of  freedom.   The  $'s  can  be  chosen  so  that  y  +  Z  3..  k.  has  an 

ij   J 

0(h   )  error  in  its  approximation  to  y(tn  +  ot.h).   In  that  case, 

s   (s) 
the  errors  in  the  estimates  of  the  h  y    will  be  proportional  to 

terms  in  h    y      +  0(h   )  and  not  involve  separate  elementary 

differentials  of  order  p  +  1.) 

It  remains  to  show  that  the  coefficients  can  be  chosen  so  that 

the  root  condition  on  Q(y)  is  satisfied,  that  is,  so  that  Q(y)  has  zeros 

in  the  positive  half  plane.   It  was  shown  in  Gear  [7]  that  Q(y)  = 

det(By  -  I)  where  B  is  the  matrix  of  beta  coefficients.   If  the  betas 

are  completely  determined  by  the  alphas  (as  happens  if  the  suggestion 

from  the  last  paragraph  is  followed),  then  Q(y)  is  completely  determined 

by  the  alphas.*   If  the  alphas  are  chosen  as  Gaussian  points  corresponding 

to  the  Butcher  2p-th  order  implicit  RK,  then  it  is  known  that  Q(y)  is  of 

degree  p  with  all  roots  in  the  positive  half  plane  (Ehle  [4]). 


*  The  zeros  of  Q  are  the  inverses  of  the  eigenvalues  of 


diag[a1?  a0,...,  a  ]A  diag[l,  1/2,  1/3,...,  l/p]A 
12       p 


-1 


where  A  = 


P"1 
1  a,   a!f 

P"1 
1  a2  a£ 


,  p-1 

1  a   a 

P   P 
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5.   PROGRAMMING  CONSIDERATIONS 

Both  the  implicit  and  explicit  techniques  can  be  used  either  to 
generate  the  Nordsieck  vector,  a  set  of  function  values,  or  derivative 
values.   However,  it  apparently  lacks  an  error  estimator,  so  what  stepsize 
should  be  used  for  the  RK  step,  and  what  stepsize  for  subsequent  multistep 
methods?  We  propose  to  copy  the  technique  followed  in  many  of  the  better 
codes  available  today,  and  use  the  estimate  of  the  error  in  a  method  of 
one  order  lower  than  that  used!   (This  is  often  called  local  extrapolation: 
it  can  be  viewed  as  forming  an  error  estimate,  and  using  it  not  only  for 
step  control  but  to  correct  the  solution.)   The  RK  Fehlberg  methods  [5] 
are  in  this  category  if  the  higher  order  answer  is  accepted,  as  are 
predictor-corrector  methods  if  the  corrector  is  one  order  higher  than  the 
predictor  and  the  predictor-corrector  difference  is  used  to  estimate  the  error, 

Since  an  estimate  of  h  y    is  available  in  the  p-th  order  method, 
this  can  be  used  to  select  the  stepsize  for  the  first  multistep  step,  and  to 
decide  if  the  RK  starter  stepsize  is  reasonable.   Assume,  for  the  moment, 
that  we  have  a  rough  approximation  to  a  step  size.   An  RK  starter  step  can 
be  taken  and  h  y    calculated.   Suppose  e  is  an  error  control  parameter. 
Three  cases  can  occur. 

(i)   ||hP  y(p)||  «  E 
(ii)   ||hP  y(p)||  »  £ 

(iii)   Neither  much  too  big  nor  much  too  small. 
In  the  first  case,  the  estimates  of  h  y    will  be  badly  tarnished  with 
roundoff  errors,  so  it  is  not  possible  to  scale.   Therefore,  the  RK  starter 
should  be  re-executed  with  a  step  about  [e/ | |hP  y  P  | |]    larger.   In  the 
second  case,  there  is  a  danger  of  the  estimates  being  badly  tarnished  by 
truncation  error,  so  the  RK  starter  should  be  repeated  with  a  smaller  h. 
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In  the  case  of  explicit  methods,  the  ratio  [e/ | |h  y    |]    may  be  about 
right;  at  least  it  should  serve  to  bring  h  y    within  the  large  tolerance 
allowed.   In  the  case  of  implicit  methods  for  stiff  equations,  an  excessively 
large  step  will  cause   |h  y    |   to  tend  to  l/|a  |,  so  it  may  be  difficult 
to  determine  the  amount  of  step  reduction.   However,  if  implicit  methods 
are  used,  some  estimate  of  the  Jacobian  must  be  available,  and  this  can  be 
used  to  get  a  crude  upper  bound  for  eigenvalues  (e.g.,  the  max  norm),  so 
that  a  step  that  is  certainly  small  enough  to  avoid  stiffness  problems 
altogether  can  be  tried  the  second  time  around  if  the  first  attempt  fails. 
In  the  third  case,  the  initial  stepsize  for  the  multistep  method  can  be 
calculated  from 

h°  "  Mc||hP/P»||' 

pl  I 

where  C   is  the  error  coefficient  of  the  multistep  method  to  be  used.   Then 
P 

the  appropriate  starting  values  can  be  calculated  with  Taylor's  series. 

What  should  be  the  criteria  for  "too  large"  and  "too  small"  in  the 

acceptance  tests.   "Too  small"  must  be  based  on  roundoff  considerations.   As 

s   (s) 
long  as  the  roundoff  errors  in  the  estimates  of  h  y    cannot  contribute  an 

error  to  the  multistep  starting  values  of  size  greater  than  £  (or  something 

a  little  smaller  for  safety),  the  RK  starter  can  be  accepted. 

"Too  big"  is  a  harder  question  to  answer.   Assuming  that 

|y|   =  1  (or  that  relative  errors  are  used),  my  intuition  is  to  require 

that   |hP  y    |   not  exceed  eP   P    .   This  is  based  on  the  idea  that  if 

(s)  -  iS   ,     ,  .       ,        ,  p+1   (p+1)  ~        ,  , 

y    =  A  ,  then  this  bound  causes  h    y      =  £,  so  that  the  errors  in  the 

RK  starter  are  of  about  the  right  size. 
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6.   DETECTING  DISCONTINUTIIES 

If  an  automatic  step  control  integrator  is  given  no  information 
about  discontinuities,  the  step  control  mechanism  can  be  very  inefficient 
in  the  neighborhood  of  a  discontinuity.   This  occurs  because  an  attempted 
step  that  straddles  a  discontinuity  yields  a  very  large  error  estimate, 
usually  causing  the  code  to  reject  the  step  and  reduce  the  stepsize 
sharply  so  that  the  step  not  only  no  longer  straddles  the  discontinuity 
but  often  so  that  it  falls  far  short  of  the  discontinuity.   Consequently, 
several  small  steps  are  taken  until  the  step  is  once  again  increased  to 
a  reasonable  size  so  that  an  attempted  step  again  straddles  the 
discontinuity.   This  process  is  repeated  with  a  large  cost  in  function 
evaluations.   Codes  have  used  a  number  of  techniques  to  reduce  the  problem. 
DIFSUB  [6],  an  early  code,  reduces  the  order  to  one  and  restarts  when  this 
type  of  difficulty  is  first  encountered,  and  then  quits  if  it  occurs  three 
times  in  succession — a  simple  but  not  too  useful  solution.   Some  codes 
resort  to  binary  division  of  the  interval  when  the  error  tests  are 
inconsistent.   (The  error  has  a  known  asymptotic  behavior  if  no 
discontinuities  are  present;  it  is  possible  to  check  for  gross  deviation 
from  such  behavior  as  an  indication  of  discontinuities.)   Binary  division 
may  be  as  good  as  any  strategy  if  no  further  information  about  the 
discontinuity  is  available.   However,  if  we  know  exactly  when  each 
discontinuity  occurs,  it  is  far  more  efficient  to  integrate  exactly  to  the 
discontinuity  and  then  restart. 

Discontinuities  can  arise  from  two  sources;   the  driving 
(t-dependent)  terms  and  the  dependent  variables.   An  example  of  a  driving 
term  discontinuity  is  smog  simulation.   The  sun's  energy  input  is  a  time- 
dependent  term  that  has  two  discontinuities  each  24  hours  as  the  sun  rises 
and  sets.   Dependent  variable  discontinuities  frequently  arise  in  simulation. 
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For  example,  a  relay  leads  to  equations  of  the  form 

if  |lR|  >  I  critical  then  switch  **■  on 

if  |lR|  <  I  min      then  switch  ■*■  off 

if  switch  =  on  then  V  =  0 
else  1=0 
where  IR  is  the  current  in  the  relay  coil,  and  V  and  I  are  the  voltages  across 
and  the  current  through  the  relay  contracts.   This  means  that  as  IR  (which 
is  a  function  of  the  dependent  variables)  increases  to  the  point  it  passes 
I  critical,  the  relay  turn  on  current,  the  coefficients  in  the  differential 
equations  change  discontinuously.   A  similar  problem  arises  in  some 
mathematical  models  that,  although  infinitely  dif f erentiable ,  exhibit  such 
sudden  changes  they  are  indistinguishable  from  discontinuities.   For 
example,  a  diode  is  sometimes  modeled  by  an  exponential  function  such  as 

t     40V     IN 

I  =  c(e    -  1) 
where  c  is  the  reverse  leakage  current.   The  constant  40,  or  some  large 
number,  causes  a  behavior  similar  to  that  shown  in  Figure  1. 


-*■  V 


Figure  1.   Analytic  function  which  behaves  like  discontinuity. 
In  such  cases  either  the  model  should  be  changed  to  a  discontinuous  model 
(e.g.,  a  piecewise  linear  function)  or  detection  schemes  which  indicate 
when  the  nature  of  the  model  changes  radically  should  be  added. 

The  basic  problem  is  to  detect  when  a  discontinuity  occurs  in 
time  so  that  when  a  step  is  about  to  straddle  it,  the  stepsize  can  be 
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reduced.   This  is  trivial  for  time-dependent  discontinuities,  so  let  us 

consider  the  problem  of  functions  of  dependent  variables.   Let  us  suppose 

that  a  discontinuity  occurs  when  an  expression  e(y)  crosses  zero.   Since 

y  is  a  function  t,  we  can  evaluate  e(y(t))  at  each  integration  step.   We 

could  use  a  set  of  past  values  e   .,  i  =  0,...,  k  from  a  multistep  method 

n-i 

to  estimate  the  time  at  which  the  next  zero  crossing  will  occur  and 
reduce  the  step  of  that  time  is  within  the  next  interval.   However,  that 
is  time-consuming  because  of  an  inverse  interpolation  in  each  interval 
and  error  prone  because  extrapolation  gives  larger  errors  than  interpolation. 
Therefore,  we  would  like  to  be  able  to  take  the  next  step  and  then  interpolate, 
This  is  possible  only  if  there  is  no  discontinuity  in  the  step.   The 
technique  employed  by  many  people  (see  Carver  [3])  is  to  delay  changing 
the  "switch"  until  the  completion  of  a  step.   Thus,  if  e(y)  is  an 
expression  such  that  a  switch  is  to  change  when  e(y)  changes  from  negative 
to  positive,  e(y)  is  evaluated  only  at  the  end  of  each  step.   If  e(y) 
remains  negative,  the  integration  proceeds  normally.   However,  if  the  value 
of  e(y)  is  found  to  be  positive  at  the  end  of  a  step,  it  indicates  a 
discontinuity  somewhere  in  that  step.   Inverse  interpolation  can  then  be 
employed  to  determine  the  t  value  at  which  e(y)  is  zero.   Interpolation 
can  then  be  employed  to  calculate  the  dependent  variables,  although  it  is 
probably  better  to  integrate  from  the  last  value  of  t,  and  possibly  iterate 
this  procedure  one  time. 
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APPENDIX 

A.   Non-Existence  of  5-Stage  Method  of  Fourth  Order 

For  a  5-stage,  fourth  order  method  to  be  possible  we  must 
be  able  to  solve 


(A.l) 


where 


A  = 


1 

2 
ai 

0 


a 


22  ai 


a. 


Ar  = 


1 

0 

0 

0 

2 

0 

0 

1 

0 

0 

0 

6 

0 

0 

3 

0 

0 

2 

0 

0 

1 

a. 


(B32  ax  +  B33  a  ) 


a. 


a 


a, 


(342  al  +  343  a2  +  344  a3} 


a, 


aA  A34 

\l  al  +   B«  4  +   B44  a3 

B«  A32  +  \k   A33 
and  T  is  a  4  x  3  matrix. 

Clearly,  the  columns  of  V   are  independent  because  the  columns  of 
the  right-han1  side  of  eq.  (A.l)  are  independent.   Let  B  be  the  4x4  matrix 


0 

a2   A32 

a3  A33 

0 

ai  A32 

B32  ai  +  333  a2 

0 

0 

333  A32 
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consisting  of  the  last  four  rows  of  A.   From  (A.l)  the  first  two  columns 

of  T  are  independent  eigenvectors  of  B  for  the  eigenvalue  zero.   Therefore, 

B  has  rank  no  greater  than  2.   If  a  =  0,  the  method  degenerates  to  a 

4-stage  method,  so  we  can  assume  a  ^  0.   Hence  the  first  row  of  B  is 

independent  of  the  other  three  because  of  the  zeros  in  the  remaining 

places  in  column  1.   Hence  the  last  three  rows  must  have  rank  at  most  1. 

Since  A   =  0,  this  is  possible  only  if  A   =  A   =  0.   Since  A   =  a  A 

72  52     62  62     1   32 

and  a  ^  0,  A   =  0,  so  that  3??  =  0.   Hence  A   =  0.   Again,  for  the 

last  three  rows  to  have  rank  1,  A_„  =  A.„  =  0.   Note  that  A_.  = 

53    63  74 

3,~  A„0  +  3,,  A„_  =  3,/  A„„  because  we  have  shown  that  A„„  =  0.   Clearly, 
43   32    44   33    44   33  32 

A  .  f   0  (otherwise  last  row  of  A  is  zero  and  it  is  impossible  to  satisfy 
eq.  (A.l)).   Hence  A  ~  ^  0.   Therefore,  since  0  =  A  ~  =  a  A   ,  a~  =  0. 
The  A  matrix  is  now 


A  = 


a 


ex. 


0    0 


0 


a 


33 

34 

1 

3 

a2 

0 

3 

a4 

0 

0 

0 

A54 

0 

0 

0 

A64 

0 

0 

0 

K, 

74 


where  A_.,  k, .    and  A_,  are  nonzero  (so  that  eq.  (A.l)  can  be  solved). 
54'   64      74 

Therefore  V        =   T        =0.   Consideration  of  the  first,  second  and  fourth 
'   41    42 

rows  of  (A.l)  as  equations  for  T        and  T        i  =  1,  2,  3  leads  to  a 
singular  system  without  a  solution,  so  a  5-stage,  fourth  order  method 
does  not  exist. 
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B.   A  9-Parameter,  6-Stage,  Fourth  Order  Method 
Again,  we  must  solve 


(B.l) 


Ar  = 


1 

0 

0 

0 

2 

0 

0 

1 

0 

0 

0 

6 

0 

0 

3 

0 

0 

2 

0 

0 

1 

where 


A  = 


a 


a. 


a 


al 
0 

0 

0 


a. 


a 


a. 


a  P 
2   2 


aiP2 


a. 


a 


a. 


°3P3 


a, 


a. 


a4PA 


ar 


ar 


a5P5 


with 
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r 


p„  = 


622ai 


(B.2) 


P~  = 


P4 


P.   = 


e32  a±   +  333  a2 


842  ax  +  643  a2  +  344  a3 


352  ax  +  353  a2  +  354  a3  +  855  c*4 


Qo  = 


B22  c^  =  a±   P2 


(B.3) 


Q3  =  332  a*  +  833  a\ 


Q4  =  342  ai  +  343  a2  +  344  a3 


1^5 


852  a1  +  353  a2  +  354  c*3  +  355  «4 


R2  =  0 


(B.4) 


R3  -  833  P2 


R4  =  e43  P2  +  344  P3 


R5  =  S3  P2  +  354  P3  +  355  P4 


The  first  six  parameters  we  will  pick  are  a.,  i  =  1,...,  5,  and  P~.   For 
the  moment  consider  P_  as  another  parameter.   These  should  be  chosen  so 
that  the  first  2x2  and  3x3  minors  of  A  are  nonsingular.   (Some 
additional  nonsingularity  conditions  will  be  imposed  later.)   Note  that 
these  determine  3?9  trivially  from  eq.  (B.2.1).   We  first  want  to  show 
that  these  determine  $   ~   and  &   ~   anc*  hence  Q  and  R_  from  (B.3. 2)  and 
(B.4. 2).   We  will  also  show  that  they  allow  R.  to  be  determined  in  terms 
of  P.  and  Q.,  i  =  4,  5.   As  before,  let  B  be  the  last  four  rows  of  A.   It 
has  rank  no  greater  than  3  because  it  is  a  4  x  5  matrix  with  two  vectors, 


F   and  T  ,  such  that 


Br.  =  o 
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T 
Therefore,  there  exists  a  row  vector  c  =  (c,  ,  c0,  c„ ,  c.)  that 

—      12    3   4 

(B.5)  c T  B  =  0 

From  the  first  column  of  B  we  have  immediately  that  c  =  0.   From  the 
second  column  we  find  that  a  c  +  a  c  =  0.   Without  loss  of  generality, 
take  c9  =  a1 ,  c„  =  -a_.   From  the  third  column  we  obtain 

(B.6)  a3  P3  c2  +  Q3  c3  +  R3  c4  =  0 

T 
If  r  is  the  third  column  of  T,   we  have  from  (B.l)  that  Bl\  =  [6,  3,  2,  1]  . 

Post  multiplying  (B.5)  by  this  column  vector  we  have 

cTBr3  =  cT[6,  3,  2,  1]T  =  0 

or 

3c_  +  2c»  +  c.  =  0  . 

2     3    4 

Thus,  c,=  -2c  -  3c9  =  2a  -  3a  .   Inserting  this  in  (B.6)  we  get 

a2  Q3  -(2a2  -  Sa^  =  a±   a3  P3 

If  we  substitute  (B.3.2)  and  (B.4.2)  in  this  we  get 

(B.7)       a^  a2  332  +  (a2  -  (2a2  -  Sa^P^e^  =  ^   a3  P3 

This,  together  with  (B. 2. 2),  forms  a  pair  of  equations  for  3^9  and  3„«. 

Consideration  of  the  fourth  and  fifth  elements  of  eq.  (B.5)  yields 

the  relation 

a.  P.  c„  +  Q.  c0  +  R.  c.  =0 
ix2    i   3    l  4 

Using  the  values  found  for  c .  we  get 

3 

ou  Q.  -  ol  a,   P, 

(B.8)  R.  =-^-i k-^-1 

x  2a  -  3a 

Note  that  the  last  row  of  A  is  now  a  linear  combination  of  the  fifth  and 
sixth  rows  consistent  with  the  right-hand-side  of  (B.l).   Hence  we  can 
now  ignore  the  last  row  of  A. 


30 

We  now  break  A  into  a  3  x  5  matrix  C  of  its  first  three  rows, 
and  the  fourth  through  sixth  rows,  D.   Then,  eq.  (B.l)  can  be  written, 
in  part,  as 


10   0 

(b.9)  cr  =   0   2   0 

_0   1   0_ 

We  have  assumed  that  the  first  3x3  principal  minor  of  A  is  nonsingular 
so  the  solution  of  (B.9)  can  be  expressed  as 


(B.10) 


r  = 


All 

A22 

0 

A22 

A22 

0 

A33 

A32 

0 

0 

0 

0 

0 

0 

0 

yil 

yi2_ 

y21 

y22 

y31 

u32 

\l 

0 

0 

y52_ 

where  jj.  =  {\i.  .},  i  =  1,  2  are  null  right  vectors  of  C  and  E  is  any  2x3 
matrix.   Note  that  since  the  first  three  column  of  C  are 


ai    a2        a2 


a,   a2   a2 


o   p2  p3 


the  values  of  A.,  are  rational  polynomials  in  P„,  linear  in  the  numerator 
and  denominator. 

The  remaining  steps  consist  of  ensuring  that  we  can  choose  E  in 
(B.10)  so  that 
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(B.ll) 


Dr  = 


0 

0 

6~ 

0 

0 

3 

0 

0 

2_ 

thus  solving  (B.l). 

Let  A.  =  {A.  .}  i  =  1,  2  be  column  vectors  and  let  e..  be  the 
-i      l  13 

elements  of  E.   From  (B.ll)  we  observe  that 

DAl  +  eH  °%  +  e21  %  =  ° 

D-2  +  612  D%  +  e22  ^2  =   ° 
Hence,  D\     and  DA  must  be  in  the  space  spanned  by  (Dy  ,  Dy  ).  Since  the 
third  column  of  (B.ll)  tells  us  that 

e13  PUj  +  e23  DM2  =  [6,  3,  2]T 

we  know  that  [6,  3,  2]  is  in  this  same  space.   If  DA  and  DA  are  independent, 

T 
which  we  now  assume,  then  they  must  also  span  this  space.   Hence  [6,  3,  2] 

is  in  span  [DA  ,  DA  ] ,  or 


(B.12) 


det 


DA1,  DA2,    3 


=  0 


This  condition  is  used  to  determine  P_  in  terms  of  the  six  parameters  a. 

3  i 

and  P?.   (In  fact,  a,  and  a^  do  not  appear  in  the  value.) 

If  we  examine  the  structure  of  jj  and  y_  ,  we  see  that  they  are 
linear  in  P,  and  P,.,  respectively,  but  otherwise  depend  only  on  the  parameters 
ot.  and  P?.   Furthermore,  y . ,  and  y   depend  only  on  the  parameters,  and  not 
on  P,  and  P-.   Because  of  this  Dy   is  linearly  dependent  on  P,  and  Q, ,  while 
Dy  is  linearly  dependent  on  P,  and  Q  .   From  the  earlier  argument  that  DA 
and  D_A„  are  in  the  space  spanned  by  Djj  and  Dy_9  we  see  that  Dy  and  Dy  must 
be  in  the  space  spanned  by  DA.  and  T)X    .      Hence 
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(B.13)  det  [DA  ,  DA2,  Dy.  ]  =0 

and 

(B.14)  det  [D,^,  DA2,  D]J2]  =  0 

Since  DA .  have  been  determined  in  terms  of  the  parameters  a.  and  P„, 
i  l      2 

(B*9+i)  i  =  4,  5,  is  a  linear  relation  between  P.  and  Q..   Thus,  we  can 
select  P,  and  P-  as  parameters,  determining  Q,  and  Q^,  and  hence  R,  and 
R,.  from  (B.8).   (In  some  cases,  vanishing  coefficients  fix  P,  and/or  P_, 
leaving  Q,  and/or  Q-  as  free  parameters.)   Now  we  observe  that  eqs.  (B.i.3) 
i  =  2,  3,  4  are  a  set  of  three  nonsingular  equations  which  determine  3,. 
j  =  2,  3,  4  uniquely.   Finally,  note  that  eqs.  (B.i.4)  i  =  2,  3,  4  are  a  set 
of  three  equations  for  $t-.  j  =  2,  3,  4,  5.   3c,-  can  be  chosen  as  a  parameter, 
determining  the  remaining  betas.   Now  the  matrix  A  is  known,  so  the 
determination  of  T   from  (B.l)  is  straightforward. 

Each  of  the  steps  given  above  is  linear  with  one  exception- 
solution  of  eq.  (B.12)  for  P~.   We  will  now  show  that  this  will  always 
have  a  real  solution,  and  thus  demonstrate  the  existence  of  solutions 
almost  everywhere  in  terms  of  the  parameters  a.,  P  ,  P  ,  P  ,  and  3,-c- 
(The  existence  of  one  solution  indicates  that  the  linear  equations  are 
nonsingular  almost  everywhere.    Such  a  solution  is  given  in  the  paper.) 

There  seems  to  be  no  simple  approach  to  demonstrating  that  (B.12) 
in  fact  reduces  to  a  quadratic  that  can  be  factored  over  the  reals  except 
to  plod  through  part  of  its  derivation.   We  will  do  this  for  the  benefit 
of  the  determined. 

If  (B.7)  and  (B.2.2)  are  solved  for  332  and  3^3,  and  the  result 


is  used  to  calculate  Q-.  from  (B.3.2)  we  get 
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Q3  =  ot1 


a2a3(ai~a2)  +  (2a2_3CXl)P2 


a2(ai"a2)  +  (2a2"3ai)P2 


Computing  the  X's  in  (B.10)  we  find  that 


h~  [ 


2       2      2 
a3  P2  -  a2  P3  ai  P3 

A      '   A 


2  „ 

'ai  P2        T 
-^  ,  0,  0]1 


(B.15) 


A2  = 


a2a3(a2-a3) 
+  2a2P3-2a3P2 


aia3^3_ai^   ala2^al~a2^ 


"  2aiP3 


+  2aiP2 


,  0,  0 


where  A  =  P_  a  a  (a  -a  )  -  P~  a  a  (a  -a  ) .   From  now  on  we  will  work  with 
A_A  and  AA?  since  this  scaling  has  no  effect  on  eq.  (B.12)  as  long  as 
A  ±   0.   Note  that  from  (B.15)  AA  and  AA2  are  dependent  (parallel)  if 

P„  a„(a  -a  )  =  P~  a„(a_-a  ).   (This  can  be  verified  by  routine  manipulation.) 

T 
Hence,  det  [ADA  ,  ADA  ,  [6,  3,  2]  ]  =  0  has  a  root  satisfying  this 

condition — which  happens  to  be  A  =  0  so  should  be  ignored.   By  computation 

we  find  that 


ADA  = 


a^(ara3)P2  -  ^(a^^ 


a1P2P3(a2-a3) 


aia2  ^ai"a2^ (a2~a3^P2P3 
a2(ai~a2)  +  (2a2~3ai)P2 


and 
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ADA2  = 


2  2  2 

a^a^o^Cc^-a  )  +  0^(0^-0^)  +  a  (a  -a  ) 


+  2a1[a2(a1-a2)P3  -  a^a^a^P^ 


a1a2a3[P3(a1-a2)  -  P^-a^ 


+  2ax   P2P3(a3-a2) 


2             2          a2a3(a  -a  )  +  (2a  -3a  )P. 
a  a  (a  -a  )P  +  a  a  (a  -a  )P  — = 

1  J  J  L     Z         l  Z     l     l     J  a2(Va2)  +  <2V3al)P2 


2a1a2(a1-a2)(a2-a3)P2P3 

2 
a  (a  -a  )  +  (2a2-3a;[)P2 


Note  that  these  are  both  linear  in  P~.   Hence  det  [ADA  ,  ADA.  ,  [6,  3,  2]  ] 
is  quadratic  in  P~.   We  have  already  seen  that  P„a~(a  -a  )  -  P„a?(a9~a  )  =  0 
is  a  zero  of  this  that  is  of  no  interest  because  A  =  0.   Hence  there  is 
exactly  one  other  real  zero  for  P~. 
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